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Ocean  prediction  systems  rely  on  an  array  of  assumptions  to  optimize  their  data  assimilation  schemes. 
Many  of  these  remain  untested,  especially  at  smaller  scales,  because  sufficiently  dense  observations 
are  very  rare.  A  set  of  295  drifters  deployed  in  July  2012  in  the  north-eastern  Gulf  of  Mexico  provides 
a  unique  opportunity  to  test  these  systems  down  to  scales  previously  unobtainable.  In  this  study,  back¬ 
ground  error  covariance  assumptions  in  the  3DVar  assimilation  process  are  perturbed  to  understand  the 
effect  on  the  solution  relative  to  the  withheld  dense  drifter  data.  Results  show  that  the  amplitude  of  the 
background  error  covariance  is  an  important  factor  as  expected,  and  a  proposed  new  formulation  pro¬ 
vides  added  skill.  In  addition,  the  background  error  covariance  time  correlation  is  important  to  allow 
satellite  observations  to  affect  the  results  over  a  period  longer  than  one  daily  assimilation  cycle.  The 
results  show  the  new  background  error  covariance  formulations  provide  more  accurate  placement  of 
frontal  positions,  directions  of  currents  and  velocity  magnitudes.  These  conclusions  have  implications 
for  the  implementation  of  3DVar  systems  as  well  as  the  analysis  interval  of  4DVar  systems. 
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1.  Introduction 

Ocean  prediction  across  the  globe  has  made  great  advances  in 
recent  decades  (Bell  et  al.,  2009).  Its  success  depends  critically  on 
the  process  of  assimilation  that  continually  corrects  a  prior  forecast 
with  recent  observations,  a  process  utilized  in  meteorology  for  dec¬ 
ades  (Kalnay,  2003).  Recent  observations,  prior  information  (a 
background  state,  denoted  as  xb)  and  dynamical  understanding 
are  combined  to  construct  an  optimal  state  estimate  as  an  initial 
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condition  for  the  subsequent  forecast  period  (Malanotte-Rizzoli, 
1996).  To  do  so,  assumptions  must  be  made  about  the  relationship 
between  the  observations  and  amongst  the  background  state  vari¬ 
ables  and  about  the  uncertainty  in  each.  In  particular,  the  error 
covariance  of  the  background  state,  denoted  as  Pb ,  is  a  key  piece 
of  information  as  Daley  (1996)  points  out:  “The  most  important 
element  in  the  statistical  interpolation  algorithm  is  the  background 
error  covariance.  To  a  large  extent,  the  form  of  this  matrix  governs 
the  resulting  objective  analysis”.  Specification  of  appropriate 
covariances  is  difficult  as  stated  by  Talagrand  (2003)  “Construction 
of  these  error  estimates  is  the  most  challenging  and  scientifically 
important  task.”  Without  confidence  in  the  formulation,  the 
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prediction  process  becomes  suspect  as  Bennett  (2002)  points  out, 
“It  is  difficult  to  develop  covariances.  It  follows  that  the  resulting 
inverse  estimate  or  analysis  of  the  circulation  also  lack  credibility.” 
Yet  formulations  of  Pb  are  rarely  tested  due  to  insufficient  observa¬ 
tions  (Brasseur  et  al.,  2005).  That  there  is  room  for  further 
improvement  is  suggested  by  isolated  examples,  such  as  the  one 
detailed  in  Section  2. 

The  ocean  is  severely  under-sampled  in  both  space  and  time, 
hindering  advances  of  many  studies  (Derber  and  Rosati,  1989). 
The  satellite  era  revolutionized  ocean  science  in  that  respect: 
Since  1992,  the  continual  presence  of  satellite  altimeters  in 
particular  provides  a  preponderance  of  information  on  ocean 
variability  across  the  globe  (Fu,  2010).  However,  the  satellite 
altimeter  constellation  remains  inadequate  for  providing  synoptic 
observations  even  of  just  the  two-dimensional  mesoscale  field  at 
the  ocean  surface  along  with  associated  fronts  and  eddies.  Le 
Traon  and  Dibarboure  (2002)  demonstrated  in  the  Gulf  of  Mexico 
that  the  lack  of  observation  results  in  substantial  errors  in 
estimating  eddy  frontal  positions.  The  ability  to  draw  clear  con¬ 
clusions  from  prior  examinations  of  possible  specifications  for 
Pb  is  limited  by  the  lack  of  data  (Lermusiaux,  2002).  As  xb  is  often 
provided  by  a  prior  model  forecast,  Gawarkiewicz  et  al.  (2011) 
estimate  Pb  by  evaluating  a  prior  forecast  with  subsequent  obser¬ 
vations  made  northeast  of  Taiwan,  thus  relying  on  a  single  sample 
to  estimate  a  subset  of  Pb.  It  is  difficult  to  obtain  a  large  enough 
number  of  independent  forecast  events  to  gain  statistical  confi¬ 
dence.  Consequently,  the  fundamental  question  remains:  To  what 
degree  are  assimilation  assumptions  leading  to  errors  in  ocean 
state  estimation? 

Our  purpose  here  is  to  systematically  evaluate  several  key 
assumptions  of  the  ocean  data  assimilation  methodology  by  utiliz¬ 
ing  the  rich  drifter  data  set  of  the  Grand  Lagrangian  Deployment 
(GLAD).  The  Consortium  for  Advanced  Research  on  Transport  of 
Hydrocarbon  in  the  Environment  (CARTHE)  deployed  295  CODE- 
type  drifters  in  the  northeastern  Gulf  of  Mexico  on  July  20-July 
31,  2012  (Poje  et  al.,  2014;  Olascoaga  et  al.,  2013;  Coelho  et  al., 
2014).  The  unprecedented  data  density  achieved  by  this  campaign 
makes  the  assessment  possible  by  sustaining  coverage  at  high 
spatial  density  over  several  mesoscale  ocean  features  and  over 
two  months.  The  primary  focus  here  is  on  the  lower  frequency 
mesoscale  circulation  in  the  deep  water  as  this  is  the  typical 
dynamical  target  of  operational  ocean  assimilation  systems.  The 
assimilation  systems  are  designed  to  constrain  the  mesoscale, 
and  thus  the  higher  frequency  solutions  have  forecast  skill 
that  is  a  function  of  the  external  forcing  and  the  dynamical 
representation. 

We  investigate  independent  perturbations  of  several  aspects  of 
the  Pb  formulation,  evaluating  the  resulting  forecasts  against  the 
dense  GLAD  drifter  observations,  which  are  not  assimilated.  Per¬ 
turbed  features  include  the  amplitude  of  the  background  error  var¬ 
iance,  horizontal  correlation  length  scales,  flow  dependent 
variations  in  correlations  and  time  decorrelation  scales.  While 
the  system  employed  here  is  an  implementation  of  3DVar,  the  find¬ 
ings  also  impact  parameter  choices  for  4DVar  systems,  which  are 
becoming  more  popular  (Cobas-Garcia  et  al.,  2012;  Janekovic 
et  al.,  2013;  Ngodock  and  Carrier,  2014).  The  detailed  study  pre¬ 
sented  here  permits  the  assessment  of  the  relative  importance  of 
each  of  the  tested  pieces  of  the  specification  of  Pb.  It  also  provides 
guidance  for  appropriate  parameter  choices,  in  particular  the  dec¬ 
orrelation  time  scales. 

After  presenting  an  example  to  motivate  the  search  for 
improved  data  assimilation  in  Section  2  and  a  synopsis  of  the  GLAD 
experiment  in  Section  3,  we  provide  details  of  the  model  setup  and 
the  experiments  with  perturbations  on  the  assimilation  back¬ 
ground  errors  in  Section  4.  The  results  are  examined  in  Section  5, 
with  discussion  in  Section  6  and  conclusions  in  Section  7. 


2.  Are  our  assumptions  suspect? 

An  example  from  the  GLAD  planning  phase  illustrates  the  short¬ 
comings  of  present  assimilation.  Daily  oceanic  condition  forecasts 
are  provided  by  numerical  model  systems  based  on  the  Hybrid 
Coordinate  Ocean  Model  (HYCOM)  and  the  Navy  Coastal  Ocean 
Model  (NCOM),  both  using  the  same  data  through  3DVar  assimila¬ 
tion  within  the  Navy  Coupled  Ocean  Data  Assimilation  (NCODA) 
system  (Barron  et  al.,  2006,  2007;  Cummings  et  al.,  2009;  Martin 
et  al.,  2009;  Rowley,  2010;  Rowley  et  al.,  2010;  Metzger  et  al., 
2010;  Smith  et  al.,  2011). 

Satellite-observed  chlorophyll  provides  an  indication  of 
Lagrangian  material  transport.  Lagrangian  Coherent  Structures 
(LCS),  which  outline  material  transport  patterns  (Haller  and 
Yuan,  2000;  Haller  and  Beron-Vera,  2012),  are  computed  from 
HYCOM  and  NCOM  model  surface  currents  and  compared  to  the 
chlorophyll  observations  (Fig.  1).  In  Fig.  la,  the  chlorophyll  plume 
from  the  high  productivity  area  around  29°N  88°W  has  spread  to 
the  southeast.  At  26°N  86.5°W,  the  plume  turns  and  extends  north¬ 
eastward,  implying  an  intense  cyclonic  feature  at  roughly  27°N 
86°W.  The  LCS  computed  from  both  models  cut  across  the  chloro¬ 
phyll  plume  at  27.5°N  87°W,  more  than  100  km  north  of  the 
observed  plume  turning.  Although  chlorophyll  is  not  an  ideal  tra¬ 
cer,  a  misalignment  of  this  magnitude  indicates  poor  agreement 
between  model  currents  and  observed  material  transport.  In  par¬ 
ticular,  neither  system  captures  the  cyclonic  turning  at  26°N. 

The  two  forecast  systems  based  on  HYCOM  and  NCOM  share  the 
same  input  data  and  data  assimilation  methodology.  The  sea  sur¬ 
face  height  anomaly  (SSHA)  along  altimeter  ground  tracks  during 
this  time  (Fig.  2)  indicates  a  cyclonic  circulation  at  27°N  86°W  that 
is  intruding  into  the  Loop  Current  Eddy  (LCE)  to  the  southwest.  An 
interpolation  of  this  data  constructed  by  AVISO  (Pascual  et  al., 
2006)  is  used  to  calculate  geostrophic  currents,  and  the  LCS  based 
on  the  geostrophic  currents  aligns  with  the  chlorophyll  imagery 
(Fig.  2). 

Given  the  3  km  model  resolution,  a  second  order  finite  differ¬ 
ence  can  reasonably  represent  first  order  derivative  wavelengths 
of  24  km  and  larger,  and  it  can  represent  nonlinear  terms  such  as 
advection  of  momentum  at  48  km  and  larger  wavelengths.  These 
are  scales  smaller  than  those  resolved  by  the  satellite  data.  The 
numerical  models  will  produce  realistic  dynamical  processes  that 
are  unconstrained  and  hence  could  exhibit  substantial  differences 
relative  to  observations.  However,  the  mislocation  relative  to  the 
chlorophyll  in  Fig.  1  is  on  the  order  of  100  km,  the  same  discrep¬ 
ancy  occurs  in  both  dynamical  systems  and  the  feature  is  resolved 
in  the  satellite  derived  LCS  in  Fig.  2.  The  conclusion  is  that  the 
altimeter  data  contain  the  essential  mesoscale  features  but  the 
data  assimilation  used  to  correct  both  models  is  faulty.  Clear  visi¬ 
ble  satellite  images  such  as  this  are  relatively  rare.  Evaluations 
are  qualitative  and  provide  only  one  snapshot.  Thus  it  is  difficult 
to  conduct  a  considered  evaluation  of  possible  error  sources  with 
only  this  data.  Fortunately,  the  drifter  observations  from  GLAD 
prove  to  be  quite  valuable. 


3.  The  GLAD  experiment 

The  drifters  are  similar  to  the  CODE  drifter  design  (Davis,  1985; 
Ohlmann  et  al.,  2001 ),  which  intends  to  measure  the  upper  1  m  aver¬ 
aged  currents.  Table  1  summarizes  the  deployment  locations,  dates, 
number  of  drifters  in  the  groups  (LSS,  SI,  S2,  Tl,  LI,  L2)  and  initial 
ocean  conditions.  Fig.  3  shows  the  number  of  active  drifters  during 
the  experiment.  Initial  group  deployment  positions  are  noted  in 
Fig.  4.  Initially,  the  reduction  in  the  number  of  drifters  is  due  to  fish¬ 
ermen  recovering  some,  and  Hurricane  Isaac  inflicts  damage  on  the 
observation  system  in  late  August.  The  slow  degradation  over  time  is 
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Fig.  1.  Satellite-inferred  chlorophyll  from  the  MODIS  Aqua  sensor  on  July  12,  2012,  in  color.  The  red  line  qualitatively  shows  the  core  of  the  chlorophyll  distribution. 
Superimposed  are  black  lines  of  the  Lagrangian  Coherent  Structures  computed  using  the  evolving  model  surface  velocities  from  the  NCOM  R1  experiment  (left)  and  HYCOM 
(right). 


Fig.  2.  The  satellite  SSHA  from  Jason-1,  Jason-2  and  CryoSat-2  (top)  during  the 
35  days  prior  to  July  19,  2012,  indicate  a  cyclonic  feature  at  27°N  86°W.  The  July  12, 
2012,  satellite-inferred  chlorophyll  image  from  Fig.  1  with  LCS  computed  from  the 
geostrophic  velocities  of  the  AVISO  composite  of  these  data  (bottom).  The  red  line  is 
the  same  as  in  Fig.  1,  which  is  used  to  identify  the  chlorophyll  core. 

due  to  the  limited  battery  life,  and  more  than  100  units  continued  to 
report  at  the  end  of  September  2012.  Drifters  report  positions  at  5- 
min  intervals  determined  by  GPS  so  that  inferred  velocities  are  quite 
accurate.  Data  gaps  occur  mostly  due  to  adverse  weather  conditions, 
which  require  careful  data  processing  to  produce  accurate  trajectory 
and  velocity  estimates.  Details  on  the  drifter  experiment  and 
processing  are  provided  in  Olascoaga  et  al.  (2013)  and  Coelho  et  al. 
(2014). 


Fig.  4  schematically  depicts  the  main  circulation  features  which 
drifters  encounter  during  August  and  September.  As  seen  in  Fig.  5, 
the  drifter  trajectories  and  model  circulation  indicate  substantial 
mesoscale  activity  in  this  region  of  the  Gulf.  The  main  features 
are  reproduced  in  the  model.  Typical  discrepancies  are  misplace¬ 
ment  of  features  such  as  August  22  where  the  cyclone  at  87.5°W 
27°N  is  observed  further  south  than  estimated  by  the  model. 
Reproducing  the  shingle  cyclones  on  the  edge  of  the  Loop  Current 
Eddy  is  particularly  difficult. 

The  drifter-inferred  velocity  spectral  content  is  summarized  by 
computing  rotary  spectra  that  contain  the  cyclonic  (counterclock¬ 
wise)  and  anticyclonic  (clockwise)  components  (Fig.  6).  The  spectra 
are  computed  by  averaging  over  all  drifters  the  amplitude  squared 
for  each  rotary  component,  and  the  amplitude  spectrum  is  the 
square  root  of  the  averaged  value.  The  spectra  for  both  rotary  com¬ 
ponents  are  dominated  by  low  frequency  motions  due  to  the  meso¬ 
scale  circulation,  both  cyclonic  and  anticyclonic  with  greater 
energy  in  the  cyclonic  counterclockwise  circulation.  Also  of  signif¬ 
icance  is  the  clockwise  energy  surrounding  one  cycle  per  day  (cpd) 
due  to  inertial  oscillations.  The  inertial  oscillations  are  a  dominant 
feature  in  the  observations  driven  by  wind  events,  and  thus  are  not 
stationary.  Wind  events  generate  inertial  oscillations  localized  in 
time  resulting  in  energy  surrounding  the  inertial  period  spread 
over  a  wide  band.  A  small  peak  near  the  M2  tidal  frequency  just 
under  2  cpd  is  also  apparent.  Both  the  inertial  and  semidiurnal 
bands  contain  high  coherence  between  model  and  drifter  spectra 
as  will  be  shown  later.  The  assimilation  systems  alter  the  energy 
at  low  frequencies  during  each  analysis  cycle.  The  variability  below 
about  0.75  cpd  is  the  focus  of  attention  to  understand  how 
assumptions  in  the  assimilation  system  affect  predictability. 

The  remainder  of  this  section  provides  a  short  description  of  the 
mesoscale  features  affecting  the  drifters  over  the  two  months.  The 
experiments  are  compared  relative  to  these  features  to  understand 
with  which  features  errors  are  associated.  The  Loop  Current  during 
July  extends  far  northwest  into  the  Gulf  of  Mexico,  and  the  LCE 
detached  just  prior  to  the  deployment,  dominating  the  area  south 
of  the  deployments.  The  LCE  northern  boundary  is  at  about  27°N 
throughout  the  experiment.  The  eastern  LCE  boundary  moves  from 
86°E  to  88°E  during  August  through  September. 

The  general  features  are  shown  schematically  in  Fig.  4.  The  First 
Cyclone  (FC)  advects  to  the  east  along  the  periphery  of  the  LCE,  and 
by  August  22,  the  FC  begins  to  merge  with  the  second  cyclone  (SC 
initially  at  28°N  87°W).  The  SC  advects  the  FC  southward,  substan- 
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Table  1 

The  summary  of  drifter  deployment  groups  indicating  the  central  location  of  the  group,  deployment  date  and  number  of  drifters  deployed. 


Drifter  group 

Initial  location 

Deployment  date  (2012) 

Number  of  drifters 

Pertinent  features  and  initial  ocean  conditions 

LSS 

DeSoto  Canyon  Area 

July  20 

20 

Large  scale  survey  covering  the  DeSoto  Canyon 

SI 

28.8°N  88.1  °W 

July  22 

90 

Initially  low  southwestward  flow 

S2 

29.2°N  87.6°W 

July  26 

90 

Spanning  observed  Mississippi  River  fresh  water  front 

Tl 

29.0°N  87.5°W 

July  29 

30 

Head  of  DeSoto  Canyon 

LI,  L2 

27.8°N  89.2°W 

July  30-31 

60 

Within  first  cyclone  (FC) 

Number  of  drifters  over  time 


Fig.  3.  The  number  of  actively  reporting  GLAD  surface  drifters  each  day  during  the 
experiment.  The  initial  deployment  events  are  noted  as  well  as  significant  events 
that  resulted  in  drifter  loss.  Slow  degradation  over  time  is  the  result  of  limited 
battery  life. 

tially  deforming  the  LCE  front,  and  the  FC  and  SC  eventually  merge. 
Northeast  of  the  SC,  an  anticyclone  (AC)  is  situated  over  the  Florida 
shelf  within  the  area  from  the  500  m  to  100  m  isobaths.  During 
August,  between  the  SC  and  AC  along  the  1000  m  isobath,  the  flow 
to  the  northwest  turns  along  the  edge  of  the  DeSoto  Canyon  to  the 
southwest.  Progressing  through  September,  this  northwestward 
flow  moves  into  deeper  waters  impinging  on  the  DeSoto  Canyon 
and  turning  to  join  the  westward  coastal  flow. 

Development  of  the  drifter  sampling  of  the  features  is  shown  in 
Fig.  5.  The  SI  drifters  advect  slowly  to  the  southwest  initially  and 
then  flow  to  the  east  between  the  FC  and  LCE.  The  drifters  become 
widely  dispersed  in  September  flowing  back  to  the  northwest.  S2 
initially  straddles  a  separatrix  with  the  western  half  of  the  drifters 
in  the  Mississippi  River  freshwater  outflow  (observed  by  the 
deploying  ship)  moving  toward  the  south  and  southwest  through 
the  DeSoto  Canyon  and  then  to  the  east  along  the  LCE.  These  drift¬ 
ers  meet  a  bifurcation  point  at  26.7°N  85.7°W  (noted  in  Fig.  4)  with 
the  majority  turning  northwest  and  the  remainder  turning  south¬ 
east.  The  eastern  half  of  S2  rapidly  flows  southeastward  into  the 
SC.  This  flow  brings  the  S2  drifters  very  close  to  the  LCE  front  on 
August  1.  On  August  8,  the  eastern  S2  drifters  move  to  the 
northeastern  side  of  the  SC.  By  August  15,  the  AC  entrains  the  east¬ 
ern  S2  drifters,  which  then  move  closer  to  the  Florida  coast,  and 
these  persist  in  the  northeastern  region  throughout  September. 
Most  of  the  LI  and  L2  drifters  are  initially  entrained  in  the  FC.  Once 
the  FC  and  SC  merge,  the  LI  and  L2  drifters  cover  a  large  area  in  the 
southeastern  quadrant  of  the  area  in  Fig.  5.  The  LI  and  L2  drifters 
gradually  then  flow  back  to  the  northwest.  The  T1  drifters  move 
southwest  through  the  DeSoto  Canyon  and  entrain  along  the  front 
between  the  merged  FC,  SC  and  the  LCE. 

4.  Assimilation  experiments 

4A.  General  system  setup  for  experiments 

The  assimilation  experiments  are  based  on  the  forecast  system 
using  NCOM  with  the  3DVar  in  NCODA,  which  are  the  operational 


Fig.  4.  General  circulation  features  (blue  lines)  during  August  2012  (top)  and 
September  2012  (bottom)  with  initial  positions  of  the  intensive  GLAD  drifter 
deployment  locations  (SI,  S2,  Tl,  LI  and  L2)  and  initial  drift  directions  (red  lines). 
The  bifurcation  point  is  first  reached  by  SI  on  August  31,  2012.  During  September 
the  drifters  are  widely  distributed. 


systems  run  at  the  Naval  Oceanographic  Office  for  high  resolution 
areas  nested  in  the  global  system.  The  domain  is  the  entire  Gulf  of 
Mexico  at  3  km  horizontal  resolution  and  50  vertical  levels.  Thirty- 
four  sigma  coordinates  are  used  above  550  m  depth,  and  sixteen  Z 
level  coordinates  are  used  below.  The  sigma  coordinate  structure 
has  higher  resolution  near  the  surface  with  the  surface  layer  hav¬ 
ing  0.5  m  thickness.  Atmospheric  forcing  is  taken  from  the  Coupled 
Ocean  Atmosphere  Mesoscale  Prediction  System  (COAMPS) 
(Hodur,  1997),  and  boundary  conditions  are  from  the  global  NCOM 
(Barron  et  al.,  2004,  2006).  Typically,  altimeter  SSHA  from  Jason-1, 
Jason-2,  and  CryoSat-2  arrive  with  24-  to  48-h  latency,  the  differ¬ 
ence  between  observation  time  and  the  assimilation  time.  The 
experiments  here  are  run  in  hindcast,  which  is  different  from 
experiments  run  in  real  time  during  GLAD  in  that  the  data  latency 
is  not  an  issue. 
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Fig.  5.  Distribution  of  GLAD  drifters  through  August  and  September.  The  drifter  tails  indicate  the  positions  over  the  prior  2  days.  The  background  color  is  the  model  sea 
surface  salinity,  and  the  black  vectors  are  the  model  surface  velocity.  The  model  results  are  from  experiment  B6,  which  has  the  lowest  errors  in  magnitude  of  velocity 
differences.  Colorbar  ranges  are  in  Practical  Salinity  Units. 


The  satellite  altimeter  SSHA  is  the  dominant  information  source 
for  constraining  the  mesoscale  field.  Altimeter-observed  SSHA  is 
used  with  the  MODAS  vertical  covariance  information  (Fox  et  al., 
2002)  to  construct  a  synthetic  temperature  and  salinity  profile, 
and  the  synthetic  profile  is  used  in  the  3DVar  assimilation.  All 
available  in  situ  data  are  also  used  in  the  experiments,  though 


there  is  typically  very  little  information  providing  synoptic  meso¬ 
scale  structure  from  in  situ  data. 

The  system  runs  a  daily  cycle  of  assimilation  and  forecast.  All 
observations  within  the  data  time  window  Tobs  prior  to  the  00Z 
analysis  time  are  analyzed  through  the  3DVar  to  produce  an  anal¬ 
ysis  increment,  which  is  added  to  the  model  system  during  an 
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Drifter  CCW  Rotary  Amplitude 


Fig.  6.  Counterclockwise  (top)  and  clockwise  (bottom)  rotary  amplitude  spectra. 
The  rotary  amplitude  squared  spectra  are  computed  independently  for  each  drifter, 
averaged  over  all  drifters,  and  the  square  root  is  plotted.  The  vertical  line  across  the 
two  plots  is  the  M2  (1.93  cpd)  frequency.  The  grey  vertical  grey  bar  (0.91-0.98  cpd) 
is  the  inertial  oscillation  range  between  27  and  29°N. 


increment  insertion  interval  Tins.  This  interval  is  6  h  for  all  experi¬ 
ments  except  for  B2,  which  uses  24  h.  The  system  divides  the 
3DVar  analysis  increment  by  the  number  of  time  steps  in  the  inser¬ 
tion  interval,  and  at  every  model  time  step  the  fractional  increment 
is  added  to  the  state.  Note  that  the  data  time  window  Tobs  and 
incremental  insertion  interval  Tins  represent  the  time  correlation 
of  the  assimilation  system.  Adding  the  analysis  to  the  initial  condi¬ 
tion  represents  a  time  correlation  function  that  is  a  Dirac  delta  in 
time,  which  is  physically  unrealistic.  Previous  experiments  show 
that  changing  only  the  initial  condition  results  in  transient  waves 
and  oscillations  that  can  persist  for  several  days  in  the  ocean, 
and  these  transients  do  not  appear  when  the  6  h  incremental  inser¬ 
tion  interval  is  used.  The  incremental  insertion  interval  represents 
a  boxcar  correlation  over  time  with  a  value  of  1  /Tins  over  Tins  and  0 
elsewhere.  Extending  the  data  time  window  Tobs  beyond  1  day 
results  in  observations  correcting  errors  with  time  scales  longer 
than  24  h. 

4.2.  Particular  assimilation  experiments 

The  assimilation  experiments  individually  perturb  different 
assumptions  within  the  background  error  covariance  function. 
The  3DVar  estimates  the  analysis  increment  Sx,  which  is  an  opti¬ 
mal  correction  to  a  background  xb  given  by  the  deceptively  simple 
equation  (Daley  and  Barker,  2001): 

<5x  =  P(,Ht[HP6Ht  +  R]_1  [y  -  HXf)]  (1) 

where  Pb  is  the  error  covariance  of  the  background  state  xb,  H  is  an 
operator  that  maps  the  background  state  to  observations  y,  and  R  is 
the  observation  error  covariance. 

A  history  of  research  into  Pb  exists  both  in  methods  for  specify¬ 
ing  an  analytic  functional  form  and  in  the  numerical  solution  pro¬ 
cess.  The  mathematical  foundation  of  the  convolution  integral 
(Gaspari  and  Cohn,  1999),  solution  through  an  implicit  diffusion 
operation  (Carrier  and  Ngodock,  2010)  and  generalized  extensions 
of  the  Gaussian  diffusion  operation  (Weaver  and  Mirouze,  2013) 
work  to  provide  state  covariance  relations  and  efficient  solutions 
based  on  historical  observations  or  heuristic  arguments  of  scales 


of  processes  such  as  the  Rossby  radius  of  deformation 

(Cummings,  2005). 

The  3DVar  here  uses  the  approach  of  Brandt  and  Zaslavsky 
(1997)  to  initially  separate  Pb  into  variance  Sb  and  a  correlation 
function  Cb  as  Pb  =  Sj/2Q,Sj/2.  A  decomposition  of  Cb  into  separable 
functions  is  then  made: 

Cfa(x,y,z,  t,  v:x',y',z',t',v') 

=  C :hb(x,y,  v,x',y,  V)Cvb(z,  v,z',  v')Cfib{x,y,  v,x!,y',  v>)  (2) 

The  correlation  between  two  variables  v  and  v\  at  horizontal 
locations  x,  y  and  x',  y\  and  vertical  positions  z  and  z'  is  the  product 
of  the  separable  components  of  horizontal  Chb,  vertical  Cvb  and  flow 
dependent  Cfdb  correlation  functions.  This  separation  of  variables  is 
relatively  typical  since  there  is  not  sufficient  a  priori  information  to 
provide  the  full  space-lagged  correlations  between  all  variables  as 
the  relations  change  throughout  the  globe.  Such  decomposition  is 
applied  in  general  circulation  models  (Derber  and  Rosati,  1989), 
the  Harvard  Ocean  Prediction  System  (HOPS)  (Lozano  et  al., 
1996),  MERCATOR  (Brasseur  et  al.,  2005),  the  Forecasting  Ocean 
Assimilation  Model  (FOAM)  (Martin  et  al.,  2007)  and  the  Climate 
Forecast  System  (CFS)  reanalysis  (Saha  et  al.,  2010).  Because  Pb  is 
a  strong  controller  of  the  solution,  additional  approaches  for  its 
specification  are  considered,  and  these  are  discussed  further  in 
Section  6. 

The  vertical  correlations  Cvb  are  computed  using  a  Gaussian 
function  with  length  scales  based  on  the  density  gradients  of  the 
background  profile.  Note  that  we  do  not  test  perturbations  of  the 
vertical  correlations.  The  altimeter  observations  are  used  in 
conjunction  with  the  correlations  based  on  historical  in  situ 
profiles  to  estimate  synthetic  temperature  and  salinity  profiles  that 
are  subsequently  provided  to  the  assimilation.  The  correlations 
change  bi-monthly  and  spatially.  Thus,  the  assimilation  process 
uses  synthetic  profiles  that  have  the  historically  observed  vertical 
correlation  impressed  on  them.  The  vertical  correlations  in  the 
assimilation  will  tend  to  smooth  the  synthetic  profile  vertical 
gradients.  Formally,  the  vertical  correlations  in  Pb  should  be  based 
on  the  observed  correlations  rather  than  constructing  synthetic 
profiles  prior  to  the  assimilation.  Since  the  synthetic  profiles  are 
based  on  historical  data,  we  believe  these  have  the  least  probability 
of  error,  and  we  have  no  other  proposed  superior  source  on  which 
to  test. 

Horizontal  correlations  Cbb  are  computed  using  a  second  order 
auto-regressive  function  (SOAR)  (Gaspari  and  Cohn,  1999)  with  a 
length  scale  based  on  the  Rossby  radius  of  deformation,  and  a 
geostrophic  coupling  component  is  included  to  relate  velocity 
and  geopotential  height.  Thus,  the  analysis  is  in  geostrophic 
balance.  The  flow-dependent  correlations  are  C/db  =  (1  +  Sf)e  sf, 
where,  sf  =  <S(SSH)/dh,<5(SSH)  is  the  difference  in  model  forecast 
SSH  between  two  observation  locations,  and  dh  is  the  specified 
flow  dependence  scale  factor,  which  is  expressed  in  centimeters. 
This  increases  the  correlations  between  points  that  have  little 
SSH  difference  and  decreases  correlation  across  SSH  gradients. 
The  SSH  field  is  used  in  the  flow  dependent  correlation  under  the 
assumption  that  the  flow  is  in  geostrophic  balance  and  directed 
along  pressure  surfaces  due  to  mesoscale  features.  Obviously  some 
processes  such  as  tides  are  not  aligned  with  this  assumption.  How¬ 
ever,  spatial  scales  for  tides  are  typically  well  separated  from  the 
mesoscale  in  the  deep  ocean.  Thus  it  is  not  expected  to  be  a  signif¬ 
icant  problem.  The  flow  field  itself  could  be  used  in  the  formula¬ 
tion,  though  some  influence  of  the  tidal  signal  would  still  contribute. 

The  3DVar  approach  assumes  all  observations  occur  at  the  anal¬ 
ysis  time,  and  thus  time  is  not  explicitly  stated  in  Eq.  (1).  Time 
scales  in  the  ocean  are  long,  and  attempts  to  take  the  time  correla¬ 
tion  into  account  mainly  consist  of  including  observations  over  a 
data  window  Tobs  that  is  long.  MERCATOR  assimilates  altimeter 
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data  covering  the  prior  7  days  in  an  assimilation  cycle  that  occurs 
every  7  days,  FOAM  has  a  daily  assimilation  cycle  with  data  used 
over  multiple  cycles  and  an  error  variance  increasing  linearly  with 
data  age,  the  Bluelink  Ocean  Data  Assimilation  System  (BODAS) 
(Oke  et  al.,  2008)  uses  an  11 -day  Tobs  to  assimilate  observations 
and  CFS  conducts  a  6  h  assimilation  cycle  using  data  in  the  prior 
10  days  with  a  weighting  based  on  data  age.  Cummings  et  al. 
(2009)  use  a  1-day  Tobs  and  daily  assimilation  cycle.  Formally, 
observations  should  be  assimilated  only  once.  Otherwise  observa¬ 
tion  error  variances  are  not  accurate.  A  long  time  between  assimi¬ 
lation  cycles  implies  that  new  data  could  improve  the  forecast  but 
is  not  having  an  impact.  While  a  one  day  assimilation  cycle  does 
use  newly  acquired  observations,  it  does  not  allow  a  long  time  per¬ 
iod  correlation.  Thus,  consideration  of  the  observation  time  win¬ 
dow  also  is  an  open  question  to  be  addressed. 

To  determine  the  impact  of  the  background  error  covariance 
components,  one  experiment  is  run  for  each  postulated  error 
source  with  changes  in  the  system  relevant  for  testing  the  specific 
sensitivity.  The  independent  changes  allow  the  results  to  be  inter¬ 
preted  more  easily.  The  main  parameter  settings  are  summarized 
in  Table  2. 

For  the  first  set  of  experiments  (denoted  with  R  for  the  first  set 
of  runs),  details  in  the  error  variance  amplitude  formulation,  hori¬ 
zontal  correlation  scale,  the  flow  dependence,  the  data  time  win¬ 
dow  and  the  incremental  insertion  interval  are  as  follows: 

R1 :  This  run  is  set  up  identically  to  the  experiment  run  in  real 
time  during  the  GLAD  deployment,  though  data  latency  is  not 
an  issue.  All  subsequent  experiments  in  the  first  set  are  devia¬ 
tions  from  this  initial  experiment.  The  data  time  window  Tobs 
includes  all  observations  obtained  within  the  prior  24  h,  and 
the  analysis  increment  insertion  interval  Tins  is  6  h.  Horizontal 
correlation  scales  are  on  average  21.2  km.  The  background  error 
variance  is  a  function  of  spatial  position  including  horizontal 
and  vertical.  It  is  computed  based  on  changes  in  the  24  h  fore¬ 
cast  state  over  the  prior  10  analysis  cycles.  Given  a  24  h  forecast 
x{  from  prior  cycle  i,  we  calculate  the  weighted  average  of 
squared  differences  Sil-iowi(x{  -  x{_i)  which  is  the  weighted 
sum  of  squared  differences  the  24  h  forecast  ocean  state  from 
different  forecast  cycles  (i  and  i-1),  where  the  weights  w, 
are  given  by  Cummings  and  Smedstad  (2013).  Areas  in  which 
observations  produce  large  differences  between  two  forecasts 
result  in  larger  variances.  This  is  an  advantage  to  the  algorithm 
as  errors  are  based  on  observation  impact  to  the  forecast.  In 
areas  that  are  not  observed  for  some  time,  the  main  contributor 
to  the  background  error  variance  results  from  model  time-var¬ 
iability.  This  is  a  drawback  due  to  the  sparse  sampling  of  the 


ocean  in  which  areas  unobserved  for  some  time  result  in  small 
differences  between  background  and  analysis.  Thus,  an  under¬ 
estimation  of  error  variance  is  expected. 

R2:  If  the  background  variance  amplitude  Sb  is  too  small,  the 
analysis  may  favor  the  background  over  the  observations,  so 
we  test  a  second  error  variance  estimation  approach  in  R2.  At 
the  conclusion  of  each  analysis,  the  forecast  error  variance  S j 
is  computed,  and  this  becomes  the  background  error  variance 
for  the  subsequent  cycle.  The  forecast  error  variance  in  R2  is 
computed  through  four  contributions: 

Sf  =  (i  -  (HrPtH  +  r)  “' )  S„  +  /  ^  w,-  (x[  -  x[_, ) 2 

k  /  i= — 10 

+  y“^wf(<5xa)i2  +  f(Sc-S())  (3) 

i=- 10 

During  the  analysis  cycle,  we  compute  the  formal  error  variance 
reduction  given  by  ^1  -  (HrPbH  +  ^jSb ,  which  accounts  for 

the  observation  distribution  and  observation  errors,  and  this 
provides  the  error  variance  at  analysis  time.  To  represent  error 
variance  at  the  subsequent  analysis  time,  this  variance  must 
increase.  Three  sources  of  error  growth  are  sequentially  added 
to  produce  the  forecast  variance  amplitude. 

The  first  error  growth  source  is  the  same  as  in  R1  and  is  added 
with  a  timescale  /  =  1/2.  The  second  growth  estimate  is  given 
by  yaX/fL-i0wi(^xi)2*  where  <5x,  is  the  analysis  increment  pro¬ 
vided  by  (1)  from  the  prior  assimilation  cycle  i,  and  w,  are 
weights  equal  to  those  in  the  first  term.  This  accounts  for  the 
recent  history  of  forecast  error  measured  by  observations.  This 
estimate  is  added  to  the  forecast  error  with  ya  =  1  /10.  Note,  this 
estimate  is  similar  to  the  estimate  in  R1  in  that  it  can  only  be 
non-zero  where  the  ocean  is  observed,  and  this  contribution 
to  the  algorithm  only  affects  areas  with  recent  observations. 
The  third  error  growth  is  given  by  yc(Sc  -  Sb),  which  compares 
the  climatological  variability  Sc  stored  in  the  Generalized  Digi¬ 
tal  Environmental  Model  (GDEM)  database  (Carnes  et  al., 
2010)  with  Sb,  and  we  increase  the  forecast  error  toward  the  cli¬ 
matological  estimate  with  yc  =  1  /20.  This  value  generally  pro¬ 
vides  an  upper  limit  on  variance  amplitude  in  the  thermocline 
in  the  absence  of  observations.  In  practice,  the  final  forecast  var¬ 
iance  amplitude  shows  the  desired  behaviors:  the  forecast  error 
estimate  is  reduced  where  observations  have  corrected  the 
model  state,  in  areas  with  observations  it  is  consistent  with 
background  minus  observation  differences,  and  it  increases 
toward  a  climatological  estimate  as  the  ocean  is  not  observed 
for  days  to  weeks.  This  forecast  error  variance  amplitude  is 


Table  2 

The  parameters  for  the  model  experiments  are  noted.  The  background  error  variance  formulations  used  in  the  experiments  are  denoted  as  R1  or  R2,  as  described  in  the  text.  A 
smaller  flow  dependent  scale  factor  results  in  smaller  horizontal  correlation  scales  perpendicular  to  the  flow  and  longer  ones  parallel  to  the  flow.  The  data  time  window  is  the 
period  relative  to  the  analysis  time  over  which  observations  are  included  in  the  analysis.  The  increment  insertion  interval  is  the  time  period  during  a  hindcast  over  which  the 
analysis  increment  is  added  to  the  system. 


Experiment 

Background 
error  variance 

Average  horizontal 
correlation  scale  (km) 

Flow  dependent  scale 
factor  dh  (cm) 

Data  time  window  Tobs 

Increment  insertion 
interval  Tins  (h) 

R1 

R1 

21.2 

12 

-24  to  0  h 

6 

R2 

R2 

21.2 

12 

-24  to  0  h 

6 

R3 

2  *  R1 

21.2 

12 

-24  to  0  h 

6 

R4 

R1 

12.6 

12 

-24  to  0  h 

6 

R5 

R1 

21.2 

12 

-7  to  0  days 

6 

B1 

R2 

21.2 

6 

-7  to  0  days 

6 

B2 

R2 

21.2 

6 

-7  to  0  days 

24 

B3 

R2 

12.6 

6 

-7  to  0  days 

6 

B4 

R2 

38.5 

3 

-7  to  0  days 

6 

B5 

R2 

21.2 

6 

-14  to  0  days 

6 

B6 

R2 

21.2 

6 

-7  to  0  days 

6 

B7 

R2 

21.2 

6 

-7  to  +7  days 

6 
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the  background  error  variance  amplitude  for  the  subsequent 
analysis  cycle.  The  parameters  of  Eq.  (3)  /,  ya  and  yc  can  be 
optimized  through  either  a  forward  perturbation  or  inverse 
estimation  approach.  Positive  results  here  indicate  increased 
skill  using  the  new  formulation,  and  thus  it  would  be  valuable 
to  determine  optimal  settings  of  the  parameters  in  the  future. 

R3:  As  a  reference  between  the  background  error  covariance  of 
R1  and  R2,  this  experiment  uses  an  error  variance  double  that  of 
Rl.  This  background  error  variance  is  larger  than  R1  though  still 
much  smaller  than  R2.  The  results  of  this  experiment  do  not 
produce  significant  improvements. 

R4:  The  LCS  results  (Fig.  1)  indicate  that  the  spatial  scales  of 
model  features  are  relatively  large  compared  to  those  in  the 
satellite  observations.  Reducing  the  horizontal  correlation  scale 
to  12.6  km  is  proposed  for  the  assimilation  system  to  represent 
the  smaller  scales  appearing  in  the  satellite  altimeter  data.  The 
results  indicate  this  did  not  provide  a  substantial  effect. 

R5:  This  experimental  case  lengthens  the  data  time  window  Tobs 
to  7  days.  In  many  prior  marine  rapid  environmental  assess¬ 
ment  exercises,  the  first  step  is  to  conduct  a  precursor  large  area 
survey  to  provide  a  large  domain  snapshot  of  dominant  regional 
features  used  to  initialize  forecasts  (Leslie  et  al.,  2008; 
Gawarkiewicz  et  al.,  2011;  Lermusiaux,  2002).  These  surveys 
are  repeated  on  a  regular  basis  throughout  the  full  experiment 
period  alternating  with  adaptive  surveys  targeting  the  local  fea¬ 
tures  of  interest.  Operators  required  accurate  local  features 
determined  by  a  full  range  of  multiple  dynamical  scales.  Assim¬ 
ilating  only  small  isolated  information  is  not  sufficient  for  fore¬ 
casts  from  local  model  nests  to  be  dynamically  balanced  with 
the  surrounding  larger  scale  features  (Robinson,  1997).  Thus, 
broad  information  is  also  required.  On  a  daily  basis,  very  iso¬ 
lated  lines  are  observed  by  satellite  in  the  Gulf  of  Mexico.  Ini¬ 
tialization  of  the  Gulf  of  Mexico  model  using  a  wide  area  ship 
synoptic  survey  is  not  feasible.  However,  during  a  7  day  period, 
satellite  altimetry  provides  wide  coverage  of  the  Gulf  and  could 
be  used  to  produce  a  fuller  domain  picture.  Given  these  consid¬ 
erations,  instead  of  using  the  24  h  of  data  as  in  Rl,  7  days  of  data 
could  provide  more  complete  information  on  the  larger  scale 
quasi-synoptic  structure.  Since  the  observation  error  levels  are 
not  changed  from  Rl,  the  repeated  use  has  the  detrimental 
result  of  effectively  increasing  data  weight.  This  limitation 
may  lead  to  some  distortions  of  the  local  dynamical  features 
and  is  addressed  later  in  experiment  B6. 

Comparison  of  the  LCS  from  the  first  set  of  experiments  shows  a 
better  qualitative  match  with  the  chlorophyll  patterns  in  R5 
(Fig.  7).  To  compute  the  LCS  we  adopt  Haller  and  Beron-Vera’s 
(2012)  proposed  geodesic  theory.  This  theory  leads  to  (attracting) 
LCS  at  any  time  t0  as  locally  minimally  stretching  strainlines 
obtained  from  a  backward  flow  computation  from  t0  to 
t0-T,T  =  30.  Such  strainlines  are  (material)  curves  tangent  to 
the  contractional  eigenvector  field  of  the  Cauchy-Green  strain  ten¬ 
sor,  given  by  C[°“r(x0)  =  [dF[°“r(x0)]  dF[°“r(x0).  Here  x0  is  the  posi¬ 
tion  of  a  fluid  particle  at  time  t0  and  F[°“r(x0)  is  the  deformation 
tensor  at  the  particle  position  at  time  t0  -  T,  which  follows  from 
solving  the  fluid  particle  motion  equation  x  =  v(x ,  t),  where 
v(x ,  t)  is  the  fluid  velocity  field,  and  the  asterisk  represents  a  con¬ 
jugate  transpose. 

The  analysis  of  this  first  set  of  experiments  is  discussed  in  Sec¬ 
tion  5.  The  results  indicate  that  the  background  error  variance  of 
experiment  R2  provides  improvement.  Comparison  to  chlorophyll 
suggests  an  improvement  in  direction  of  currents  in  experiment 
R5,  and  the  implication  is  that  there  is  a  long  time  correlation  in 


R5 


Fig.  7.  LCS  structures  from  the  experiment  R5  superimposed  on  MODIS  chlorophyll 
image  on  July  12,  2012.  The  change  in  the  data  time  window  Tobs  in  R5  produces  the 
best  agreement  from  the  first  set  of  experiments.  The  red  line  is  the  same  as  in 
Fig.  1,  which  is  used  to  identify  the  chlorophyll  core. 


the  background  errors  that  should  be  addressed.  However,  drifters 
suggest  experiment  R5  contains  shortcomings  in  the  magnitude  of 
the  currents.  Based  on  this,  a  second  set  of  experiments  is  defined 
denoted  with  a  B  as  the  second  set  of  runs; 

B1 :  This  is  a  combination  of  the  modified  background  error  var¬ 
iance  used  in  R2  and  the  lengthened  data  time  window  of  R5 
using  the  7  day  data  time  window.  All  subsequent  experiments 
in  this  second  set  are  a  perturbation  from  this  reference.  The 
flow  dependent  scale  factor  dh  is  decreased  in  the  B  experi¬ 
ments,  which  results  in  increasing  observation  correlation  in 
the  direction  of  constant  geopotential  and  decreasing  correla¬ 
tion  in  the  perpendicular  direction. 

B2:  This  is  the  same  as  Bl,  except  that  the  increment  insertion 
interval  Tins  is  increased  from  6  h  to  24  h.  Using  a  24  h  interval 
would  be  expected  to  produce  more  realistic  correction  to  the 
background  state,  though  a  shorter  increment  insertion  interval 
requires  less  hindcast  time.  Thus,  there  is  motivation  to  reduce 
the  increment  insertion  interval  from  the  computational  cost 
perspective.  Results  show  this  parameter  has  a  negligible  effect. 
B3;  The  horizontal  correlation  scale  is  again  tested  in  combination 
with  the  long  time  data  window  by  reducing  the  value  to  1 2.6  km. 
The  results  again  did  not  show  significant  improvement. 

B4:  The  horizontal  correlation  scale  in  the  direction  of  constant 
surface  geopotential  is  increased  in  this  experiment  by  chang¬ 
ing  the  flow  dependent  scale  factor  dh  and  changing  the  hori¬ 
zontal  correlation  scale.  The  flow  dependent  scale  factor 
increases  horizontal  correlation  scales  in  the  direction  of  the 
constant  geopotential  and  decreases  scales  in  the  direction  per¬ 
pendicular.  At  the  same  time,  the  horizontal  correlation  scale  is 
increased  to  38.5  km.  The  result  of  these  two  effects  is  to  main¬ 
tain  the  same  horizontal  correlation  scale  perpendicular  to  the 
flow  and  greatly  extend  the  horizontal  correlation  scale  parallel 
to  the  flow.  The  purpose  is  to  extend  observation  influence  fur¬ 
ther  in  the  direction  of  flow.  This  experiment  did  not  result  in  a 
positive  improvement. 

B5;  This  experiment  doubles  the  data  time  window  to  14  days 
from  the  7  days  used  in  Bl.  The  results  do  not  demonstrate 
increased  skill. 

B6;  Experiment  Bl  uses  observations  multiple  times  because  of 
the  7-day  data  time  window,  and  this  process  is  formally  incor- 
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Fig.  8.  Background  error  standard  deviation  Sb  of  190  m  temperature  on  2012  August  21  through  25  (top  to  bottom  rows),  from  experiments  (left)  R1  and  (right)  R2. 
The  circled  area  indicated  on  August  23  contains  reduced  variance  from  observations  on  August  22  (Fig.  10).  Note  the  different  colorbar  ranges  for  R1  (0.25-0.75  °C)  and  R2 
(0-2  °C). 
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rect.  Given  the  observation  error  levels,  using  observations 
more  than  once  results  in  the  solution  matching  the  observa¬ 
tions  more  than  it  should,  and  the  analysis  is  more  susceptible 
to  noise  in  the  data.  A  correct  method  to  use  observations  once 
and  provide  a  long  time  correlation  would  be  to  use  a  data  time 
window  of  1  day  and  extend  the  incremental  insertion  interval 
over  multiple  days.  As  previously  discussed,  the  system  divides 
the  analysis  increment  by  the  number  of  time  steps  in  the  incre¬ 
ment  insertion  interval.  Thus  the  difference  between  inserting 
over  one  day  and  inserting  over  7  days  is  dividing  the  analysis 
increment  by  a  factor  of  7.  This  experiment  is  the  same  as  B1 
(daily  cycling  with  a  7  day  data  time  window  and  6  h  incremen¬ 
tal  insertion  interval)  with  the  exception  that  the  analysis  incre¬ 
ment  is  divided  by  the  number  of  days  in  the  data  time  window. 
This  is  intended  to  emulate  using  observations  in  one  analysis 
with  an  analysis  increment  insertion  interval  of  7  days  with  a 
correlation  that  is  constant  in  time  throughout  the  7  days.  This 
is  not  an  entirely  correct  emulation  since  the  analysis  incre¬ 
ment  influence  of  single  observations  changes  from  one  day 
to  another  as  the  background,  and  thus  observation  minus 
background  changes. 

B7:  As  a  reanalysis,  it  is  possible  for  data  in  the  future  to  affect 
the  analysis  time  as  well  as  data  in  the  past.  This  experiment 
uses  a  data  time  window  that  is  -7  days  to  +7  days  relative  to 
the  analysis  time.  Results  did  not  indicate  a  significant  benefit. 


5.  Results 


5.1.  Background  error  variance  amplitude  and  time  correlation 


The  background  error  variance  amplitude  Sb  formulation  is  a 
strong  controller  of  the  solution.  Fig.  8  provides  a  comparison  of 
Sb  over  several  days  computed  by  the  error  variance  formulations 
in  R1  and  R2.  The  amplitudes  are  roughly  a  factor  5  larger  in  R2 
compared  to  R1  (note  different  colorbar  ranges  in  Fig.  8  for  R1 
and  R2).  The  estimated  error  variance  amplitude  Sb  is  evaluated 
relative  to  the  observed  errors  in  the  background  field  (b-o).  At 
each  depth  d,  with  Nd  subsurface  observations  during  the  experi¬ 
ment,  an  estimate  of  background  error  variance  in  temperature 
Sd  is  computed  as 


Sd 


(4) 


where  the  background  field  is  interpolated  to  the  observation  loca¬ 
tion.  For  comparison,  the  estimated  background  error  variance  Sb  in 
experiments  R1  and  R2  are  sampled  at  the  observation  locations 
and  averaged  over  the  experiment  time.  Fig.  9  compares  the  square 
root  of  the  error  variance  from  R1  and  R2  to  the  estimate  Sd.  The 
variance  values  of  Sb  should  be  less  than  those  of  S^since  observa¬ 
tion  representation  errors  and  noise  contribute  to  Sd  in  addition  to 
the  background  errors.  Fig.  9  shows  the  Sb  in  experiment  R2  is  a  bet¬ 
ter  representation  than  experiment  R1  when  compared  to  observa¬ 
tion  minus  background  variance  at  most  depths.  The  typical 
thermocline  and  halocline  depths  are  about  200  and  250  m  respec¬ 
tively  in  the  eastern  Gulf  of  Mexico.  Results  from  R2  indicate  the 
background  error  variance  Sb  computed  by  the  algorithm  in  Eq. 
(3)  provides  a  good  representation  of  the  estimated  Sd,  though 
the  algorithm  provides  a  slight  over  estimate.  One  reason  for  the 
apparent  over-estimation  may  be  that  the  observations  are  mainly 
the  synthetic  profiles  generated  from  the  altimeter  data.  The  syn¬ 
thetic  profiles  will  contain  less  variance  than  real  profiles. 

The  R2  increment  fields  (Fig.  10)  reflect  the  location  of  satellite 
tracks.  In  R2,  the  background  error  variance  reduction  based  on 
observation  locations  and  error  levels  is  computed  and  subtracted 
from  Sb,  and  then  the  error  growth  terms  are  added  according  to 


(3).  This  becomes  the  background  variance  for  the  successive  analy¬ 
sis  cycle  on  the  next  day,  so  that  observation  impacts  on  Sb  are  seen 
in  the  background  variance  on  the  following  day.  One  example  is 
indicated  by  the  circles  in  Figs.  8  and  1 0.  On  August  22,  observations 
over  the  western  Gulf  of  Mexico  (Fig.  10)  produce  background  error 
variance  reduction  on  August  23  (Fig.  8).  The  background  error  var¬ 
iance  on  successive  days  (Fig.  8)  shows  the  impact  of  relaxing  toward 
climatological  variability  in  areas  that  have  no  observations  as  the 
error  variance  under  prior  observations  gradually  increases. 

The  second  important  factor  is  the  correlation  time  scale  of  the 
background  errors.  The  increments  are  compared  between  R2  and 
B1  (Fig.  10),  which  both  use  the  same  Sb.  B1  has  a  7  day  data  time 
window  vs.  1  day  for  R2.  For  example,  on  August  21  there  are 
increments  along  the  satellite  track  transiting  from  20.5°N  83°W 
to  30.5°N  86°W  (encompassed  by  the  box  in  Fig.  10).  In  experiment 
R2,  the  increments  occur  along  this  track  only  on  August  21.  In 
experiment  Bl,  the  increments  appear  on  August  21,  and  subse¬ 
quent  days  also  show  the  increments  with  amplitude  diminishing 
in  time  over  3-4  days.  If  background  errors  are  not  correlated  over 
long  times,  the  second  use  of  the  data  would  have  a  much  smaller 
increment  and  the  increment  would  be  more  noise  than  signal. 
Additionally,  the  increments  in  the  case  of  Bl  are  generally  smaller 
than  the  increments  due  to  the  same  data  in  R2,  indicating  the 
background  field  is  closer  to  the  observations  in  Bl. 

The  increment  fields  on  August  21  of  experiments  B1-B7 
(Fig.  11)  indicate  the  effects  of  the  increment  insertion  interval 
(B2),  the  flow  dependent  correlations  (B3  and  B4),  the  data  time 
window  (B5  and  B7),  and  emulating  the  long  time  correlation 
(B6).  The  main  improvement  relative  to  the  drifters  occurs  in  B6 
in  which  the  analysis  increments  are  divided  by  the  data  time  win¬ 
dow  resulting  in  much  smaller  values. 

5.2.  Evaluation  through  GLAD  observations 

The  unassimilated  GLAD  drifters  serve  as  independent  evalua¬ 
tion  data.  Each  experiment  is  sampled  at  the  latitude  and  longitude 
location  of  each  drifter.  The  velocity  vectors  are  computed  every 
3  h,  resulting  in  over  100,000  observations.  The  3-h  sampling 
observes  the  same  mesoscale  eddies  for  quite  some  time,  thus 
the  number  of  independent  observations  is  not  the  same  as  the 
number  of  observations.  The  drifters  do  cover  a  wide  range  of  dif¬ 
ferent  eddies  and  associated  fronts  during  the  deployment  time 


Fig.  9.  Square  root  of  the  specified  background  error  variance  averaged  over  all 
observations  during  the  experiment  from  R1  (solid  black)  and  R2  (dashed  black) 
compared  to  the  RMS  difference  of  observations  minus  background  from  R1  (solid 
red)  and  R2  (dashed  red). 
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Fig.  10.  Analysis  increments  on  August  21  through  25  2012  (top  to  bottom  rows)  for  (left)  R2  with  a  1  day  data  time  window  and  (right)  B1  with  a  7  day  data  time  window. 
Highlighted  areas  are  discussed  in  the  text.  These  are  the  corrections  to  the  model  state  computed  by  the  3DVar  system.  Colorbar  ranges  are  in  °C. 
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Fig.  11.  Analysis  increments  on  August  21,  2012  from  experiments  perturbing  properties  relative  to  B1  including:  increment  insertion  interval  (B2),  the  flow  dependent 
correlations  (B3  and  B4),  the  data  time  window  (B5  and  B7),  and  emulating  the  long  time  correlation  (B6).  Colorbar  ranges  are  in  °C. 


period  (Fig.  5),  thus  the  number  of  independent  observations  is 
high  during  the  60  day  experiment  period. 

An  initial  comparison  is  made  through  the  rotary  coherence  and 
phase  difference  spectra  (Fig.  12).  The  model  experiments  indicate 
high  coherence  over  the  mesoscale  band  from  0  to  0.75  cpd  in  both 
clockwise  and  counterclockwise  spectra  with  a  phase  difference 
between  ±20°.  Typically,  R2  has  higher  coherence  than  R1  in  this 
band,  R5  is  lower  than  Rl,  and  B6  is  near  the  top  coherence.  As  a 
reference,  geostrophic  currents  from  the  AVISO  sea  surface  height 
maps  are  computed,  and  the  statistics  for  these  currents  are  also 
presented  in  Fig.  12  and  subsequent  analyses. 

The  inertial  oscillations  in  the  0.9-1. 1  cpd  range  in  the 
clockwise  spectra  show  high  coherence  in  the  models,  though 
phase  differences  are  higher.  The  inertial  oscillations  are  typically 
forced  by  surface  wind  stress  events.  The  accuracy  in  the  inertial 
band  is  a  measure  of  wind  forcing  accuracy  and  dynamical  system 
response  in  the  surface  layer.  In  addition,  coherence  is  high  at  the 


semidiurnal  frequency  just  below  2.0  cpd,  and  this  is  due  to  tidal 
forcing  at  the  boundaries  and  tidal  potential  forcing  in  the  interior. 

The  high  frequency  (greater  than  0.75  cpd)  response  is  primar¬ 
ily  a  forced  response,  and  the  interest  at  hand  is  in  the  low 
frequency  response  that  is  affected  by  the  assimilation  process. 
To  focus  on  this  variability,  model  velocities  are  interpolated 
through  a  bilinear  interpolation  to  the  latitude  and  longitude  loca¬ 
tions  of  the  drifters.  Then,  all  the  time  series  of  velocities  from  the 
models  and  those  inferred  from  the  drifters  are  passed  through  a 
Butterworth  filter  with  a  0.5  cpd  cutoff.  Thus  the  subsequent  error 
characterizations  are  for  the  low  frequency  only. 

Error  histograms  over  all  data  are  computed  (Fig.  13).  These  are 
the  fraction  of  observations  within  each  bin  treating  magnitude 
and  direction  errors  independently  rather  than  binning  error  in 
the  vector  difference.  Thus,  skillful  system  error  distributions  should 
be  clustered  around  0°  direction  error  and  small  magnitude  error. 
The  distribution  for  Rl  provides  the  fraction  of  all  observations  that 
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Fig.  12a.  Counterclockwise  (top)  and  clockwise  (bottom)  rotary  coherency  spectra.  A  filter  with  frequency  width  of  0.1  cpd  is  applied.  The  colored  lines  are  the  coherency 
spectra  with  vertical  scale  on  the  right,  and  the  shaded  background  is  the  amplitude  spectrum  from  the  drifters  as  in  Fig.  6  with  the  vertical  scale  on  the  left. 
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Fig.  12b.  Counterclockwise  (top)  and  clockwise  (bottom)  rotary  phase  difference.  A  filter  with  frequency  width  of  0.1  cpd  is  applied.  The  colored  lines  are  the  phase  difference 
spectra  with  vertical  scale  on  the  right,  and  the  shaded  background  is  the  amplitude  spectrum  from  the  drifters  as  in  Fig.  6  with  the  vertical  scale  on  the  left. 


fall  in  each  (magnitude,  direction)  bin.  The  majority  of  errors  is 
clustered  around  0°  direction  and  is  less  than  0.2  m/s.  For  more 
direct  comparison,  the  error  distribution  of  R1  is  subtracted  from 
the  histograms  computed  from  other  experiments,  and  the  color 
bar  range  is  changed  to  highlight  the  difference  in  fraction  of  obser¬ 
vations  between  the  experiments  and  R1 .  Areas  of  blue  indicate 


fewer  occurrences  than  Rl,  and  areas  of  red  indicate  more  occur¬ 
rences.  An  experiment  is  performing  better  than  Rl  if  there  is  a  dis¬ 
tribution  of  blue  away  from  the  center  and  red  near  the  center  and 
near  0°  as  in  the  distributions  of  R2  and  B6. 

A  16  element  matrix  fully  describes  cross  correlations  between 
observed  and  modeled  vectors  (Crosby  et  al.,  1993),  and  a  range  of 
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Fig.  13.  Error  histograms  distributions  relative  to  GLAD  inferred  velocity  observations  from  R1  (top  row).  Skillful  system  error  distributions  should  be  clustered  around  0° 
direction  error  (toward  top  of  the  figure)  and  small  magnitude  error.  Other  experiments  are  presented  as  differences  in  distribution  minus  the  distribution  of  R1  (second  row) 
AVISO  and  R2,  (third  row)  R5  and  B6.  An  experiment  is  performing  better  than  R1  if  there  is  a  distribution  of  blue  away  from  the  center  and  red  near  the  center  and  near  0°  as 
in  the  distributions  of  R2  and  B6. 


statistics  are  proposed  using  subsets  of  these  elements.  To  help 
illustrate  the  results  succinctly,  we  compute  cumulative 
distributions  that  show  the  fraction  of  the  observations  with  errors 
less  than  a  prescribed  value.  It  is  useful  to  examine  these  across 
three  statistics:  (1)  the  difference  between  observed  and  modeled 
magnitude  of  velocity  vector  (Fig.  14),  (2)  difference  in  direction 
measured  in  degrees  (Fig.  15),  and  (3)  the  magnitude  of  the  differ¬ 
ence  in  observed  and  modeled  velocity  vector  (Fig.  16).  The  cumu¬ 
lative  probability  distributions  of  errors  are  similar,  and  small 
differences  are  important.  To  better  highlight  the  small  differences, 
the  cumulative  distribution  for  experiment  R1  is  subtracted  for 
each  of  the  statistics,  and  the  difference  between  R1  and  each 
experiment  is  also  presented  in  the  figures.  In  comparing  two 
experiments,  the  distributions  with  more  occurrences  at  low  error 


values  are  performing  better.  Because  these  are  cumulative  distri¬ 
bution  curves  as  a  function  of  error  level,  the  higher  (and  more 
positive)  the  curve  at  low  error  levels,  the  better  the  performance. 
Only  the  experiments  that  demonstrate  significant  impact  beyond 
the  initial  experiments  R1  and  B1  are  shown,  which  are  R2,  R5  and 
B6.  The  geostrophic  currents  computed  from  direct  analysis  of 
SSHA  by  AVISO  are  included  as  a  reference. 

R2  is  performing  overall  better  than  R1  with  a  higher  fraction  of 
occurrences  at  low  error  levels  across  all  statistics.  This  is  confir¬ 
mation  that  the  Sb  of  R2  is  a  better  representation  than  the  Sb  value 
of  Rl.  R3,  which  doubles  the  Sb  of  Rl,  has  little  impact  on  the 
results.  Similarly,  R4,  which  reduces  the  horizontal  correlation 
scale  using  the  same  error  variance  amplitude  as  Rl,  has  little 
impact.  R3  has  a  distribution  very  similar  to  Rl,  and  R4  performs 
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Cumulative  distribution 


Difference  in  cumulative  distribution 


Fig.  14.  (Top)  cumulative  distributions  of  the  difference  in  magnitude  of  velocity 
from  the  model  experiments  and  the  GLAD  drifter  inferred  velocities.  (Bottom)  to 
better  visualize  the  differences,  the  distribution  of  R1  is  subtracted  from  all  the 
experiments.  Distributions  with  more  occurrences  at  low  error  values  are 
performing  better.  Thus  B6  is  the  best  performing  experiment  by  this  metric. 


slightly  worse.  The  cumulative  distributions  are  not  shown  for 
these  for  brevity  and  clarity. 

R5  is  the  first  experiment  to  include  a  long  data  time  window 
(7  days).  Though  the  initial  LCS  comparison  in  Fig.  7  indicates  an 
improved  qualitative  comparison  relative  to  the  advected  chloro¬ 
phyll,  the  error  distribution  statistics  from  the  GLAD  drifters  show 
poorer  performance  than  Rl.  The  errors  in  magnitude  of  velocity 
are  the  primary  cause,  and  this  is  also  reflected  in  the  magnitude 
of  the  difference  between  observed  and  modeled  R5  velocities. 
The  positive  comparison  to  the  chlorophyll  distribution  points 
out  the  risk  involved  in  evaluations  based  on  a  small  sample  set 
or  metric,  which  represents  one  instant  in  time  and  only  an  evalu¬ 
ation  of  flow  direction  compared  to  the  chlorophyll. 

The  comparison  between  R2  and  Bl,  which  both  use  the  same 
Sb  but  have  different  data  time  windows  (Fig.  10)  implies  that 
background  errors  are  correlated  in  time.  This  should  be  expected. 
Many  estimates  of  the  error  covariance  structure  start  with  the 
assumption  that  the  scales  of  the  errors  are  on  the  order  of  the 
scales  of  the  features  themselves.  The  ensemble  optimal  interpola¬ 
tion  of  BODAS  (Oke  et  al.,  2008)  uses  ocean  model  ensembles  to 
build  the  spatial  covariance  relations.  The  11 -day  assimilation 
cycle  of  this  system  also  points  to  the  use  of  a  long  data  time  win¬ 
dow.  Since  ocean  mesoscale  features  have  time  scales  of  weeks 
(Jacobs  et  al.,  2001 ),  it  should  be  expected  that  errors  are  correlated 
in  time.  The  experiments  altering  the  increment  insertion  interval 
(B2),  flow  dependent  covariance  (B3,  B4),  and  additional  changes  in 


Cumulative  distribution 


Error  in  direction  (degrees) 


Fig.  15.  (Top)  cumulative  distributions  of  the  difference  in  direction  of  velocity 
from  the  model  experiments  and  the  GLAD  drifter  inferred  velocities.  (Bottom)  to 
better  visualize  the  differences,  the  distribution  of  Rl  is  subtracted  from  all  the 
experiments.  R2  is  performing  best  with  more  occurrences  at  low  errors.  B6  is 
performing  better  than  Rl. 


the  data  time  window  (B5,  B7)  do  not  significantly  improve  the 
results. 

The  experiment  combining  the  new  Sb  with  the  emulation  of 
long  time  correlation  is  B6,  and  this  shows  the  best  performance 
with  the  most  occurrences  of  small  velocity  magnitude  error 
(Fig.  14),  more  occurrences  of  low  errors  in  direction  than  Rl 
though  not  as  good  as  R2  (Fig.  15)  and  the  best  performance  in 
terms  of  most  occurrences  of  low  errors  in  magnitude  of  velocity 
difference  (Fig.  16).  The  use  of  a  long  time  window  forces  the 
model  solution  to  match  the  data  within  the  time  window.  R5  suf¬ 
fers  from  the  fact  that  the  solution  must  match  the  average  of  the 
observations  used  in  the  analysis.  As  the  observation  time  window 
increases,  the  analysis  weakens  horizontal  gradients  resulting  in 
correspondingly  weaker  currents.  Experiment  B6  avoids  this  prob¬ 
lem  by  applying  the  analysis  increment  over  a  long  period  with 
reduced  amplitude  that  allows  the  dynamical  system  to  maintain 
the  horizontal  gradients  with  increased  current  magnitudes. 

Fig.  17  provides  a  summary  spatial  distribution  of  the  RMS  vec¬ 
tor  difference  magnitudes  in  Rl,  that  is 

-  u„,)2  +  ( -  vjf,  where  umi  and  uoi  are  model 
and  observed  velocities  respectively.  These  statistics  are  computed 
in  1/8°  bins,  and  only  bins  with  more  than  30  samples  are  shown. 
The  data  are  separated  for  August  and  September.  Areas  of  larger 
error  are  situated  in  the  fronts  of  the  mesoscale  features,  which 
move  through  August  and  September.  The  area  along  the  northern 
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Fig.  16.  (Top)  cumulative  distributions  of  the  magnitude  of  velocity  difference 
between  the  model  experiments  and  the  GLAD  drifter  inferred  velocities.  (Bottom) 
to  better  visualize  the  differences,  the  distribution  of  R1  is  subtracted  from  all  the 
experiments.  The  B6  experiment  has  the  highest  number  of  occurrences  at  low 
errors  and  therefore  is  performing  best. 


edge  of  the  Loop  Current  Eddy  has  large  errors  associated  with  the 
First  Cyclone  (FC)  advecting  southeastward.  The  FC  can  be  seen  in 
the  drifter  tracks  of  Fig.  5  on  August  15  at  27°N  88°W.  The  FC  being 
less  than  50  km  across  can  be  marginally  resolved  in  the  models, 
though  it  is  smaller  than  what  could  be  expected  to  be  captured 
by  the  satellite  observations.  The  boundaries  of  the  AC  on  the  Flor¬ 
ida  shelf  are  apparent  as  are  the  currents  to  the  northwest  between 
87  and  88°W  27°N. 

The  spatially  distributed  errors  of  the  perturbed  covariance 
experiments  are  shown  in  Fig.  18,  and  the  circulation  features  of 
Fig.  4  are  superimposed.  For  a  clearer  comparison,  the  difference 
in  the  statistics  for  each  experiment  and  those  of  R1  are  shown 
for  August  and  September.  Blue  areas  show  errors  lower  than  R1 , 
and  red  areas  are  errors  higher  than  R1 .  Experiment  R5  indicates 
the  degradation  in  velocity  performance  in  August  with  increased 
errors  in  the  FC  along  the  northern  LCE  edge,  the  bifurcation  point 
at  85.5°W  27°N,  the  Florida  shelf  AC  and  an  overall  increase  in 
errors  in  September  in  the  northwestward  flow.  Experiments  R2 
and  B6  both  indicate  substantial  improvement  over  R1  and  R5. 
During  August,  R2  predicts  the  FC  along  the  LCE  front  with  lower 
errors  than  B6.  Relative  to  experiment  R2  in  August,  B6  improves 
the  bifurcation  circulation  at  85.5°W  27°W,  the  Florida  shelf  AC, 
and  the  SC  circulation  at  87°W  28°N.  During  September,  B6  gener¬ 
ally  decreases  overall  error  levels  throughout  the  region  relative  to 
R2.  Thus  experiment  B6  indicates  improved  performance  across 
more  mesoscale  features  than  the  other  experiments. 


5.3.  Interpretation  in  the  context  of  4DVar 

A  4DVar  assimilation  can  be  reduced  through  simplifying 
assumptions  to  produce  the  3DVar  used  here.  It  is  useful  to  exam¬ 
ine  the  assumptions  to  provide  insight  to  the  aspects  3DVar  is 
neglecting  and  how  those  assumptions  are  reflected  in  the  results. 
Consider  the  4DVar  solution  that  minimizes  errors  relative  to  a 
time-evolving  background  state  xb,  observations  y  and  the  dynam¬ 
ics  A,  i.e.,  that  minimizes  the  cost  function  J: 

J  =  I  <)xtP,7'  5x  +  I  SxtAtW'ASx  + 1  (H<5x  -  d)TR  1  (H<Sx  -  d)  (5) 

where  <5x  is  the  analysis  increment,  Pb  is  the  background  error  covari¬ 
ance,  A  is  the  tangent  linearization  of  the  dynamical  system  around 
the  background  state,  W  is  the  dynamical  error  covariance,  H  is  the 
linearized  observation  operator,  d  =  (y  -  H(xb))  is  the  misfit  of  the 
observations  to  the  background  state  xb  and  R  is  the  observation 
error  covariance.  A  particular  choice  of  covariances  in  the  4DVar 
can  reduce  it  to  the  3DVar.  Consider  the  4DVar  applied  over  the  entire 
time  period  of  interest  (60  days).  All  model  variables  at  all  locations 
and  over  all  time  steps  are  concatenated  to  form  a  single  vector 
representing  the  model  state  trajectory  xb  through  the  entire  time 
period.  Assume  that  the  temporal  structure  of  background  error 
covariance  Pb  is  taken  to  be  a  sum  of  Dirac  delta  functions  at  the  anal¬ 
ysis  times  of  the  3DVar  3(t  -  tai ),  where  tai  is  the  time  of  the  ith  3DVar 
analysis.  This  results  in  Pb  being  block  diagonal  with  blocks  providing 
the  interrelations  amongst  variables  only  at  each  analysis  time. 
Assume  that  observations  occur  only  at  the  analysis  times.  The  solu¬ 
tion  to  (5)  can  be  obtained  by  an  inversion  conducted  at  each  succes¬ 
sive  analysis  time.  Assume  the  solution  is  a  strong  constraint  except 
at  the  analysis  times,  and  thus  satisfies  the  dynamics.  Once  the  ith 
inversion  is  conducted,  the  time-evolving  state  is  computed  simply 
by  integrating  the  analysis  forward  to  the  next  analysis  time.  Under 
these  conditions,  (5)  represents  the  cost  function  for  a  cycling  3DVar. 

The  block  diagonal  form  of  Pb  carries  strong  implications.  Back¬ 
ground  errors  beyond  the  analysis  time  are  assumed  to  be  uncor¬ 
related.  This  is  convenient  for  solving  the  problem  in  that  the 
analysis  can  be  computed  for  each  cycle  independently.  However, 
this  formulation  precludes  consideration  of  errors  correlated 
beyond  one  assimilation  cycle.  Fig.  10  shows  that  error  correlations 
extend  out  to  days.  Thus,  when  considered  from  the  4DVar  per¬ 
spective,  the  cycling  interval  defined  in  a  3DVar  imposes  an 
assumption  limiting  the  background  error  time  correlation. 

A  4DVar  system  can  use  a  long  time  window  in  the  analysis,  and 
its  contribution  to  reducing  the  analysis  error  covariance  is  com¬ 
puted  theoretically  for  linear  systems  by  Lewis  et  al.  (2006).  The 
authors  also  provide  a  detailed  derivation  of  the  difference  between 
the  analysis  error  covariances  of  the  4DVar  and  the  Kalman  filter  (of 
which  3DVar  is  a  special  case),  which  includes  a  proof  that  the  Kal¬ 
man  filter  solution  error  covariance  is  always  greater  than  or  equal 
to  the  4DVar  solution  error  covariance.  The  difference  lies  in  the 
4DVar  ability  to  propagate  the  observation  information  further  in 
space  and  time.  The  implicit  restriction  imposed  by  the  3DVar 
assimilation  cycle  results  in  an  increase  in  error  levels. 

The  improved  formulation  for  the  variance  amplitude  dem¬ 
onstrated  here  can  also  be  used  in  the  4DVar.  However,  as  dis¬ 
cussed  above,  error  levels  in  the  3DVar  are  expected  to  be  higher 
than  the  4DVar.  While  a  similar  algorithm  for  estimating  the  back¬ 
ground  error  variance  may  be  used  with  a  4DVar,  the  amplitudes 
must  be  smaller. 

For  global  daily  predictions  that  are  computationally  intensive, 
a  4DVar  solution  is  not  feasible.  Thus  there  is  continued  need  to 
understand  both  the  3DVar  and  4DVar  approaches.  Analytic  for¬ 
mulations  do  not  represent  many  complex  relations  within  Pb,  thus 
Fu  (2012)  estimates  Pb  based  on  historical  model  runs  for  methods 
such  as  ensemble  optimal  interpolation  (EnOI)  and  shows  benefits 
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Fig.  17.  (b)  Number  of  observations  in  1/8°  bins  during  August  and  September  2012  and  RMS  magnitude  of  the  vector  velocity  error  of  experiment  R1  during  August  (a)  and 
September  (c).  Errors  are  shown  only  in  bins  with  more  than  30  observations.  Error  colorbars  are  in  m/s. 


over  3DVar.  However,  basing  Pb  entirely  on  model  covariances  can 
introduce  errors  due  to  biases  and  drifts  in  the  models.  A  stationary 
Pb  as  implemented  by  Oke  et  al.  (2010)  avoids  problems  in  model 
drifts,  but  it  may  still  be  influenced  by  biases.  This  motivates  the 
hybrid  approach  in  which  Pb  is  specified  as  a  linear  combination 
of  model  covariances  and  covariances  based  on  functional  formula¬ 
tions  as  in  Eq.  (2)  (Gao  et  al.,  2013). 

Another  perspective  on  the  time  correlation  comes  from  the 
Rapid  Environmental  Assessments  in  which  a  large-scale  survey 
provides  an  initial  state.  Attempting  to  correct  a  system  based  on 
under  sampled  features  results  in  poor  performance.  On  a  daily 
basis,  the  satellite  coverage  is  unable  to  resolve  mesoscale  features 
as  can  be  seen  in  Fig.  10,  which  shows  days  with  very  little  satellite 
data.  Thus  it  is  not  possible  to  accurately  correct  these  features 
using  only  daily  data.  Data  spanning  several  days  can  begin  to 
resolve  features.  However,  a  negative  consequence  is  that  a  3DVar 
analysis  results  in  a  state  averaged  over  the  data  time  window.  In 
addition,  using  the  observations  multiple  times  without  correcting 
the  observation  error  level  over-constrains  the  solution  to  the 
observations.  These  are  the  problems  causing  poor  performance 
in  R5.  The  temporal  averaging  results  in  weak  horizontal  gradients 
and  thus  weak  currents.  The  R5  distribution  of  error  in  magnitude 
(Fig.  14)  and  distribution  of  magnitude  of  velocity  difference 
(Fig.  16)  suffer  due  to  this  problem.  The  R5  distribution  of  direction 
errors  (Fig.  15)  does  not  suffer  as  badly.  There  are  roughly  equal 
occurrences  of  direction  errors  less  than  35  degrees  as  for  Rl. 
Experiment  B6  mitigates  the  deficiencies  in  using  observations 
multiple  times  by  emulating  the  effect  of  a  long  time  correlation 
through  dividing  the  analysis  increment  by  the  number  of  days 
in  the  increment  insertion  interval.  As  shown  in  Fig.  11,  the 
increments  in  B6  are  very  small  and  applied  over  many  days. 


The  results  from  B6  indicate  much  improved  performance  over 
R5  in  the  distribution  of  difference  in  magnitude  and  magnitude  of 
the  velocity  differences.  B6  has  the  highest  fraction  of  occurrences 
at  the  low  error  levels  of  all  the  experiments  (including  those  not 
shown).  The  directional  errors  in  B6  show  advancement  over  Rl 
at  low  values  but  are  somewhat  degraded  relative  to  R2. 


6.  Discussion 

The  GLAD  experiment  shows  it  is  now  possible  to  begin  testing 
underlying  assumptions  of  ocean  assimilation  methods,  which  has 
not  been  extensively  possible  before.  While  persistently  sampling 
several  mesoscale  features  at  high  resolution  over  two  months, 
the  data  set  is  still  just  one  experiment  and  exists  only  in  the  Gulf 
of  Mexico.  The  results  and  conclusions  should  be  transferable  to 
other  regions  and  to  global  mesoscale  prediction  as  well.  However, 
predictive  skill  varies  within  different  dynamical  regions  such  as 
the  Gulf  Stream  and  Kuroshio  currents  vs.  the  South  China  Sea 
(Metzger  et  al.,  2014b).  Because  the  Gulf  of  Mexico  is  enclosed,  pre¬ 
diction  skill  is  higher  in  this  basin  using  the  same  metrics  (Metzger 
et  al.,  2014a).  Also,  mesoscale  prediction  skill  is  highly  dependent 
on  the  quantity  and  coverage  of  altimeter  observations  (Ananda 
et  al.,  2006;  Smedstad  et  al.,  2003;  Jacobs  et  al.,  2014). 

The  GLAD  data  point  to  additional  underlying  considerations  in 
future  direction  of  ocean  assimilation.  Many  ocean  data  assimila¬ 
tion  methods  are  formulated  with  the  objective  to  predict  the 
ocean  mesoscale  as  developed  through  the  GODAE  program 
(Cummings  et  al.,  2009).  This  objective  is  feasible  given  the  altim¬ 
eter  observations,  which,  though  sparse  in  space  and  time,  are 
capable  of  sampling  a  portion  of  the  mesoscale  wavenumber  and 
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Fig.  18.  Differences  between  RMS  magnitudes  of  the  vector  velocity  error  computed  for  R1  and  other  experiments.  Error  levels  of  square  root  of  mean  of  vector  error 
magnitude  squared  differenced  from  experiment  R01.  Only  bins  with  more  than  30  observations  are  shown.  August  and  September  2012  are  shown  separately.  The 
circulation  features  of  Fig.  4  are  superimposed.  Colorbar  units  are  in  m/s. 
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frequency  spectrum  (Traon  and  Dibarboure,  1999).  The  long  time 
scales  of  the  ocean  result  in  a  mesoscale  feature  being  sampled 
several  times  over  a  decorrelation  time  of  one  month  or  more. 
One  area  that  connot  yet  be  constrained  with  in  assimilation  is 
the  submesoscale  which  is  not  resolved  by  the  altimeter 
observations,  has  large  Rossby  number  and  high  ageostrophic 
currents  (Capet  et  al.,  2008;  Mensa  et  alM  2013;  Zhong  and 
Bracco,  2013). 

Three  fundamental  capabilities  must  exist  to  address  the  sub¬ 
mesoscale.  The  first  is  observations.  GLAD  is  a  demonstration  that 
sufficient  coverage  and  persistence  of  observations  at  scales  below 
the  mesoscale  are  achievable.  Other  planned  observation  systems 
include  the  Surface  Water/Ocean  Topography  (SWOT)  satellite 
mission  planned  to  provide  full  global  2D  sea  surface  height  imag¬ 
ing  with  a  resolution  of  kilometers  (Durand  et  al.,  2010;  Fu  and 
Ferrari,  2008). 

The  second  capability  lies  in  multiscale  data  assimilation.  Basic 
research  has  yielded  algorithms  that  solve  the  assimilation  prob¬ 
lem  sequentially  over  scales  (Brandt  and  Zaslavsky,  1997;  Choi 
et  al.,  2008;  Haley  and  Lermusiaux,  2010;  Lermusiaux,  2002). 
Application  to  the  submesoscale  requires  a  first  multiscale  itera¬ 
tion  assimilating  the  mesoscale  information.  This  problem  has 
been  well  studied.  The  second  multiscale  iteration  should  then 
address  the  submesoscale. 

The  third  capability  must  provide  a  new  Pb  for  the  submeso¬ 
scale  that  removes  assumptions  present  in  Pb  for  the  mesoscale 
iteration.  Amongst  these  are  geostrophic  balance  relations 
between  the  mass  and  velocity  fields  that  are  inappropriate  for 
the  submesoscale  as  well  as  vertical  relations  in  the  water  column. 
Solving  the  submesoscale  iteration  of  multiscale  assimilation 
requires  revisiting  the  background  error  covariance  Pb  that  will 
have  distinctly  different  horizontal,  vertical  and  temporal  struc¬ 
tures  than  the  mesoscale.  Once  again,  we  will  be  required  to  test 
new  Pb  formulations  for  the  submesoscale. 


7.  Conclusions 

Our  analysis  shows  that  background  error  variance  amplitude 
and  time  correlation  are  the  most  sensitive  components  within 
the  background  error  covariance  Pb.  The  extraordinarily  dense 
GLAD  data  set  provides  substantially  greater  information  to  inves¬ 
tigate  the  relative  importance  of  assumptions  in  the  components 
of  Pb.  This  is  a  problem  not  previously  addressable  due  to  sparse 
observations.  The  assimilation  systems  are  used  to  emulate  the 
long  time  correlations  by  using  data  over  several  prior  days  in 
each  daily  cycle.  When  using  data  multiple  times,  the  modeled 
currents  are  weak  due  to  the  analysis  providing  a  best  estimate 
to  match  an  average  over  7  days.  This  is  alleviated  by  dividing 
the  analysis  increment  fields  by  the  data  time  window  length 
to  more  properly  weight  the  data.  The  results  pave  the  way 
toward  implementations  that  can  properly  account  for  the  time 
correlations  important  in  3DVar  and  are  more  explicitly  apparent 
and  controllable  in  4DVar  through  the  analysis  time  window. 
While  4DVar  can  take  into  account  a  longer  influence  time  of 
observations,  care  must  be  taken  to  ensure  the  analysis  window 
is  sufficiently  long  to  account  for  the  long  time  correlation  in 
the  ocean. 

A  study  such  as  this  would  not  be  possible  without  an  exten¬ 
sive  set  of  observations  covering  a  significant  number  of  meso¬ 
scale  features  for  a  sustained  time  period.  Prior  hypothesized 
formulations  for  the  background  error  variance  amplitude  and 
time  correlation  are  rejected  in  favor  of  new  understanding.  Such 
data  sets  become  reference  points  that  are  used  throughout 
future  research.  The  GLAD  experiment  points  the  direction 
toward  the  future. 
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